% ========================================================
% 主程序：使用Cholesky分解求解二维泊松方程
% 方程形式：-aΔu = f 
% ========================================================

close all
% 定义输入参数
n = 20;          % 单方向内部节点数
A_1 = 10;        % 右端项第一高斯源幅值
A_2 = 10;        % 右端项第二高斯源幅值
h = 1/(n+1);     % 步长定义

% ----------------------------------------
% 计算所有矩阵和向量
% ----------------------------------------
% 生成n*n维刚度矩阵
S = DiscretePoisson2D(n);

% 通过Cholesky分解矩阵A = L*L^T
[L] = chol(S,"lower");

%% 生成系数矩阵a(x_i,y_j)
C = ones(n,n);  % 此处为常系数1的情况

%% 计算载荷向量f（双高斯源）
f = zeros(n^2,1);
for j = 1:n
    for i = 1:n
        pos = n*(i-1)+j;
        term1 = exp(-((i*h-0.25)^2/0.02 + (j*h-0.25)^2/0.02));
        term2 = exp(-((i*h-0.75)^2/0.02 + (j*h-0.75)^2/0.02));
        f(pos) =A_1*term1 + A_2*term2;
    end
end

% ----------------------------------------
% 使用高斯消元法求解线性系统
% ----------------------------------------
% 系统方程：1/h² (C*L*L^T) u = f

% 求解线性系统
b = f; % C为全1矩阵，b = f./C(:) 简化为b = f
v = L \ b;       % 前向替换
w = L' \ v;      % 后向替换
u = h^2 * w;     % 缩放解

% ----------------------------------------
% 结果可视化
% ----------------------------------------
% 将解向量排列为网格格式（边界节点为0）
Z = zeros(n+2);
Z(2:end-1, 2:end-1) = reshape(u, n, n)';

%% 绘图
x1 = 0:h:1; y1 = 0:h:1;

figure(1)
surf(x1, y1, Z)
view(2), colorbar
xlabel('x_1'), ylabel('x_2'), zlabel('u(x_1,x_2)')
title(['解场分布图 N=', num2str(n)])

figure(2)
surf(x1, y1, Z)
colorbar
xlabel('x_1'), ylabel('x_2'), zlabel('u(x_1,x_2)')
title(['三维解场视图 N=', num2str(n)])
